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ABSTRACT 

This paper presents a framework for developing com- 
putationally unified numerical algorithms for solving non- 
linear equations that arise in modeling various problems in 
mathematical physics. The concept of computational unifi- 
cation is an attempt to encompass efficient solution proce- 
dures for computing various nonlinear phenomena that may 
occur in a given problem. For example, in Computational 
Fluid Dynamics (CFD), a unified algorithm will be one that 
allows for solutions to subsonic (elliptic), transonic (mixed 
elliptic-hyperbolic), and supersonic (hyperbolic) flows for 
both steady and unsteady problems. The objective of the 
work reported in this paper is manyfold: 1) development of 
superior unified algorithms emphasizing accuracy and effi- 
ciency aspects; 2) development of codes based on selected 
algorithms leading to validation; 3) application of mature 
codes to realistic problems; and 4) extension/application of 
CFD-based algorithms to problems in other areas of math- 
ematical physics. The ultimate objective is to achieve in- 
tegration of multidisciplinary technologies (stealth, propul- 
sion, aeroelasticity, • • •) to enhance synergism in the design 
process through computational simulation. 

The paper presents specific unified algorithms for a hi- 
erarchy of gasdynamic equations (full potential, Euler, and 
Navier-Stokes) and their applications to a wide variety of 
problems. Also included are extensions of the CFD meth- 
ods to two other areas: 1) electromagnetic scattering, and 
2) laser-material interaction accounting for melting. 

INTRODUCTION 

Along with rapid strides in algorithm and code de- 
velopment, the increasing power of super-minicomputers, 
supercomputers, and graphics workstations is rapidly ad- 
vancing the state of the art of computational simulation 
of problems in mathematical physics. One area that is set- 
ting the pace is Computational Fluid Dynamics. Other dis- 
ciplines such as electromagnetic scattering, semiconductor 
device/process modeling, material characterization, etc., 
are starting to benefit from the CFD experience. 

Modern vehicle concepts such as the Advanced Tacti- 
cal Fighter (ATF) attempt an effective compromise between 
the transonic maneuver and supersonic cruise conditions. 
Multiple design considerations of this type impose strin- 
gent constraints on the aerodynamic shape of the vehicle to 
achieve high buffet-free lift performance with reduced trim 
drag. The recent resurgence of the hypersonics program 
through the National Aerospace Plane (NASP) project also 
demands analysis and design of vehicles with requirements 
to fly through the entire Mach number range (subsonic to 
hypersonic) requiring increasingly sophisticated nonlinear 
methods to better understand various gasdynamic flow pro- 
cesses. 


The Navier-Stokes equations best represent the phys- 
ics of nonlinear flow. However, limitations in memory and 
execution speed of present-day supercomputers restrict the 
routine use of Navier-Stokes methods. For wider applica- 
tion of CFD in the aerospace industry, cost-effective meth- 
ods based on less exact forms of gasdynamic equations, such 
as the Euler and full potential equations, are still attrac- 
tive. The objective of the work reported in this paper is 
to develop, for all speed regimes, efficient, accurate, and 
robust nonlinear methods for equations ranging from the 
simple full potential to the complex Navier-Stokes. De- 
velopment of such a spectrum of hierarchical capability is 
critical for efficient and cost-effective design of aerospace 
configurations. The general philosophy of numerical design 
through progression of increasingly sophisticated nonlinear 
tools is illustrated in figure 1, and represents a summary of 
the numerical design experience at Rockwell covering the 
HiMAT, forward swept wing, SAAB, and Air Force/Navy 
Research Technology contract studies^"**. 

Referring to figure 1, in designing a configuration, lin- 
ear theory^’® is first used to establish candidate optimum 
thickness, twist, camber, and variable camber deflections 
at supersonic speeds. Second, nonlinear methods (full po- 
tential and Euler) are employed to capture embedded shock 
waves at transonic^” and supersonic^^”*^^ conditions and 
weaken the wave system through parametric redesign. 
Boundary layer analysis^^ and Navier-Stokes codes^^“^® 
are subsequently used to assess the flow quality of the non- 
linear inviscid design. The extent of separation in particu- 
lar is evaluated, and a subsequent redesign is performed to 
minimize its extent. 
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Fig. 1. Progression of improved design through nonlin- 
ear analysis. 
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At Rockwell, the computational activity is carried out 
on several fronts. 1) Algorithm development . Under this 
task, several algorithmic issues such as higher order space/ 
time accuracy, efficiency, multizone gridding concepts, 
multigrid cycling, upwind total variation diminishing 
(TVD) schemes, vectorization concepts, implicit/explicit 
methods, etc., are stressed. The primary thrust of this ac- 
tivity is computational unification of methodologies encom- 
passing efficient solution procedures for computing various 
flow phenomena occurring across the Mach number range. 
2) Code development . Under this task, selected algorithms 
(based on the study from task 1)) with potential to mature 
into a production code capability are further developed and 
undergo extensive validation. Some of the issues stressed 
in this phase are computational efficiency through code 
vectorization, user orientation/documentation, code trans- 
portability, graphics interface, and user training. 3) Appli- 
cation . At this stage, codes from task 2) that have matured 
into a production code with established user confidence, are 
applied to study a wide range of realistic problems to bet- 
ter understand flows over existing configurations (Shuttle 
Orbiter, B-lB, etc.) as well as to design aerospace config- 
urations for the next generation (ATF, NASP, Transatmo- 
spheric Vehicles (TAV), etc.). 4) Extension of CFD 
methods to non-CFD problems . Under this task, a host of 
problem areas in mathematical phy?ks that are governed 
by appropriate partial differential equations is dealt with. 
The techniques developed for studying problems in Compu- 
tational Fluid Dynamics are well suited for studying prob- 
lems in electromagnetic scattering, laser-material interac- 
tion, and semiconductor device/process modeling, just to 
name a few. 

Computational Fluid Dynamics is rapidly advancing. 
Its methods are beginning to influence how problems are 
and can be effectively solved in other disciplines. As this 
process of spreading the wealth of CFD knowledge to other 
areas continues, the CFD discipline is expected to play a 
key role in the future, together with state-of-the-art com- 
puters, in integrating multidisciplinary technologies to en- 
hance synergism in the design process through computa- 
tional simulation. 

The paper presents a brief summary of some of the 
unified algorithms developed for various gasdynamic equa- 
tions, along with their applications to many fluid dynamic 
problems. Also presented are some applications of CFD 
methods to non-fluid dynamic problems. More details on 
the algorithmic aspects of the unified concept can be found 
in the references. 

This paper represents a collection of work performed 
by many researchers in the Computational Fluid Dynamics 
Department at Rockwell International Science Center. 

EQUATIONS IN CONSERVATION FORM 

For many problems in mathematical physics, the phys- 
ical process to be modeled is governed by an appropriate 
set of linear or nonlinear partial differential equations. For 
example, many fluid dynamic processes are governed by 
the Navier-Stokes equations, the electromagnetic scatter- 
ing from objects is modeled by Maxwell’s equations^^, and 
problems in semiconductors are governed by Van Roos- 


broeck^* equations involving the nonlinear coupling be- 
tween the electrostatic potential and the electron/hole den- 
sity. 

In general, many of these equations naturally lend 
themselves to a conservation form representation given by 

Qt d" d" d" Gz 0 (I) 

where the dependent variable vector Q, and the fluxes E, 
F, and G take on different forms depending on the equation 
being modeled. The form of Q, E, F, and G for the full 
potential, Euler, Navier-Stokes, and the Maxwell equations 
are presented in the subsequent sections. Application of 
eq. (1) to many realistic problems requires a coordinate 
transformation to properly represent the physical domain 
of interest and to aid in the boundary condition treatment. 

Under the transformation of coordinates implied by 
T = <,^ = = ri(t, x,y,z),i = f(<,x,y,z), 

eq. 1 can be recast in the conservation form given by 

+ + =0 , (2a) 

where 

F= + , 

G=%Q+^E+^F + ^G . 

where, in turn, J is the Jacobian of the transformation 

./ = 9(^,t/,f)/a(x,y,z) (2c) 

and 

d- d- ^z^t) 

rit = -(Vx^T d- Tiyyr d- rjzZr) (2d) 

U ~ ~{^x^r d" d- ^z^t) 

Associating the subscripts with the direc- 

tions, a numerical approximation to eq. (2a) may be ex- 
pressed in the semi-discrete conservation law form given by 

(Qy,fc,i)r d- (Ej+if2,k,l - 

d" (^y,fc+l/2,/ ~ Ej^k~lf2,i) (^) 

-h (Gj^k,l+i /'2 ~ 1-1(2) = 0 

where F, F, G are numerical or representative fluxes at the 
bounding sides of the cell for which discrete conservation is 
considered, and Qj\k,i is the represent ^ive conserved quan- 
tity (the numerical approximation to Q) considered conve- 
niently to be the centioidal value. The half-integer sub- 
scripts denote cell sides and the integer subscripts the cell 
itself or its centroid. 
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The semi-discrete conservation law given by eq. (3) 
may be regarded as representing a finite volume discretiza- 
tion if the following associations are made: 


Qj\kj = (4a) 

where V is the volume of the cell under consideration; 


V J )r 


(46) 


y±i/3 

- 1/2,/ - 1/2), (A; + 1/2,/ - 1/2), 

(k + 1/2,/ + 1/2), (A: - 1/2,/ + 1/2)}, ±1/3 

/ flx,y,z \ _ 

V J /*±l/2 

n^.vAU - 1/2,/ - 1/2), (; - 1/2,/ -f 1/2), 

U + 1/2,/ + 1/2), (; + 1/2,/ - l/2)}fc±i/3 
/ _ 

^ J //±l/2 

n.^yAU - 1 / 2 , k - 1 / 2 ), (; + 1 / 2 , k - 1 / 2 ), 

{j + 1/2, A: -h 1/2), (j - 1/2,/: -f- l/2)}/ii/2 

(y) = “ (y )y±i/a(®r)y±i/2 

~('j))±l/‘ 2 iyr))±l /2 ~ (y)y±l/2(*r)y±l/3 

““(^)fc±l/2(yr)/t±l/2 - {^)fc±l/2(^r)fc±l/2 

-(^)i±1/2(Vt)i± 1/2 - (y)l±l/2(^r)l±l/2 

In the above, rix^y^z are the x^y^z components of the repre- 
sentative normals to the surface formed by the four points 
a,6, c, d implied in nx,y,z(cLib,c,d). These four points are 
not necessarily coplanar. Also, ( 2 :^, J/r, ^r)y±i/ 2 ,fc±i/ 2 ,f±i /2 
are the x, y, z components of the appropriate cell-face rep- 
resentative velocities. These describe the motion of the cell 
face and will be zero for a stationary grid. In the following, 
we can use the notation n* to describe the representative 
cell-face normal velocities: 


(4c) 


(««)y±i/2 - ^y^ 
(^t]k±i/2 = (y) 
(^t)l±l/2 - ( 


i±l/3 

fc±l/2 

l±lj2 


(id) 


The evaluation of the volume, cell-face normals and 
cell-face normal velocities (metrics) are presented in 
Ref. 14. 

Within this framework of a finite volume representa- 
tion, the concept of a unified algorithm/solver addresses 
two issues: 1) representation of the numerical fluxes E, F, 
and G to account for different physical phenomena to be 
encountered in the problem being modeled (for example, in 
fluid dynamics, a unified flux representation will allow for 


subsonic, transonic, and supersonic flow situations for both 
steady and unsteady, including proper transition through 
shocks and sonic rarefaction); and 2) numerical issues of 
solving eq. (3). Under the numerical issue, a unified algo- 
rithm will be one that performs both the space and time 
integration within the logic of a single solver. A unified 
solution treatment will allow one to consider a wide class 
of problem areas within the capability of a single code. For 
example, in fluid dynamics, a unified solver will perform 
space marching for supersonic flows (supersonic flow direc- 
tion is treated as time like) and allow time marching for 
subsonic, transonic, and unsteady flows. 

The objective is to solve eq, (3) for the dependent vec- 
tor Q. After incorporation of proper flux representation, 
the discrete form of eq. (3) can be written as 

R(Q)=0. (5) 


If Q is known at a known neighborhood state, denoted by 
then solution to eq. (5) can be written as 

-<?•) = - W) (6) 

where in general, is a differential operator. Many 
numeric^ algorithmic issues such as implicit, explicit, re- 
laxation, approximate factorization, algorithm unification, 
etc., come into play in the modeling of the differential op- 
erator 1^. Issues such as higher order accuracy, proper 
upwinding, etc., come into as well as in the modeling of 
the right hand side /?(<?*). For a unified code that accounts 
for both space and time marching, one option is to split the 
operator in the form 


where 


dR 

dQ 




( 7 ) 


Ln=L^[(r, 0 ,r>] , - L, [(r, ^ f] • 


Equation (7) represents a double approximate factorization 
in the (f/, f ) plane with relaxation in the ^-direction as- 
sumed to represent the predominant flow direction. The 
grouping (r, ^) in the and operator represents a col- 
lection of terms involving time and ^ derivative terms. For 
time marching, the time-step-size Ar is chosen to main- 
tain the stability and accuracy of the operator, eq. (7). For 
space marching, Ar is usually set very large and the oper- 
ator Lfi becomes L,,(^,»;) and = Lj(^,^) representing f 
as the marching direction. Space marching along ^ is pos- 
sible only if the equation is hyperbolic with respect to that 
direction. A code that is based on the unified solver will 
include the following options: 

II = Lf{i,T)L^(fl,T)L({i,T) - Triple approxi- 

mate factoriza- (8a) 
tion with time 
marching 

II = L, [(t, ^), .j)] L, [(r, 0. f] - Double approxi- 

mate factoriza- (8b) 
tion with time 
marching 
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if = I>iU,n)I>sU,v) - Space marching 

along ( setting (8c) 
Ar 00. 

If one employs an upwind differencing, can also be rep> 
resented by a Gauss- Siedel relaxation maintaining diagonal 
dominance. More discussions on these ideas can be found 
in Refs. 12, 19, and 20. 

For time marching, Q* is usually set to be Q” as a 
first guess where Q”’ is the solution at the previous time 
plane. For space marching, Q* is initially set to be Qj-i 
representing the solution at the previous space marching 
plane in the ^ direction. Starting from the initial guess 
for Q*, eq. (6) is iterated to convergence driving || Q — 
Q* II to some preset small value at every time or space 
marching plane. Usually, this process might involve only a 
few iterations. 

The issue of numerical flux representation is dealt with 
in the subsequent sections for a variety of equations, namely 
1) full potential, 2) Euler, 3) Navier-Stokes, 4) Maxwell, 
and 5) incompressible Navier-Stokes, representing the 
laser-material interaction. 

Full Potential Equation 

The full potential equation represents the inviscid, ir- 
rotational, and isentropic flow. In spite of these assump- 
tions, this form of the gasdynamic equation is widely in 
use for analyzing complex configurations at transonic and 
low supersonic Mach numbers. As long as the shocks are 
weak (Mach number normal to a shock surface less than 
1.3 to 1.5), the full potential isentropic shocks will be in 
agreement with the Rankine-Hugoniot jump conditions. 


is a differential operator. 

Exjuation (3) also requires evaluation of E, F, and G 
at various sp^ial half node points. As mentioned earlier, 
E represents E appearing in eq. (2a). 

The concept of developing a unified full potential 
scheme stems from a grope^definition for the numerical or 
representative fluxes E, F, G at cell interfaces derived from 
the theory of characteristic signal propagation. Depending 
on the type of flow at the cell interface (subsonic, transonic, 
or supersonic), the fluxes are properly defined employing an 
upwind bias to eliminate numerical or spurious (unphysical) 
oscillations by satisfying entropy conditions (no expansion 
shocks). 

Based on the characteristic system at a cell interface 
in the time-space domain, the following different flux rep- 
resentations are made. 

Subsonic at cell face j + 1/2 (U < Cy/an) 

Ej+i /2 = ^y+ 1/2 (zero biasing) 

Transonic at j + 1/2 {U < Cy/a^ , 9 > c) 

^j+i /2 = -^y+i /2 (to be defined later) (12) 
Supersonic at j -f 1/2 {p > c^an) 

Ej^ij 2 = Ej_ij 2 (upwind biased flux) 

In eq. (12), c is the speed of sound and an — + + ^ 2 )* 

The transonic flux E is defined in terms of an upwind biased 
density based on flux biasing. Define 


Referring to eq. (1), the full potential equation takes 
the form Q = E = pu, F ~ pv, and G = pu;, where p 
is the density and u, v, and w are the Cartesian velocities. 
All the quantities p, u, v, and w are expressible in terms 
of a single scalar function the velocity potential. Using 
Bernoulli’s law, the density p is given by 

= l-^AC[2^r + (C + 6)^e jgj 

and Uy V , and W are the contravariant velocities. 

Referring to eq. (3), modeling of the time term 
will require time linearization for density to express AQ in 
terms of A^ = <f> — <p* . The density linearization is given 
by 



- 1 r, , , f C/ d V d 


where Q — \/U'^ + + W^. 

The quantity [pq)~ appearing in eq. (13) is defined to 


{pq) =pq- p*q* 


\f q>q* 
ii q<q* 


The quantities p*q*, p*, and q* represent sonic values 
of the flux, density, and total velocity, respectively. These 
sonic conditions are given by (using the density and speed 
of sound relationships) 

/ - Ur - UU, - UUr, - Hth) 

w ) - 

2 -*”oo 


/>* = (<7* Moo) 


2/(7-1) 


Note that for steady flows, the sonic conditions p* and 
q* are only a function of the freestream Mach number, and 
for a given flow they are constants. For unsteady flows, p* 
and q* need to be computed everywhere due to the presence 
of 4>f and other unsteady terms in eq. (15). 
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In the final discretized form, the full potential equation 
is written in the form of eq. (6) with = (^ — ^*) 3 ls a 
single unknown at a grid point. 

Some results are presented to illustrate the unified full 
potential capability in computing subsonic, transonic, and 
supersonic steady/unsteady flows. 

Fall Potential Resnlts 

Supersonic Flows — Supersonic flows are computed using 
the marching option within the unified solver. 

Figure 2 shows the surface gridding along with cross- 
plane field grid points for a typical advanced generic fighter. 
The body-fitted grid is generated at every marching plane 
using standard elliptic grid solvers. Figure 3 shows pres- 
sure contours at two different axial stations at Moo = 1-fi, 
a = 4.94°. The crossplane geometry in figure 3 clearly 
shows the fuselage, vertical tail, wing, and the flow through 
nacelle along with the wake cut behind the trailing edge of 
the wing. Figure 4 shows pressure correlation between the 
computations and experimental data at two different span 
stations. Table 1 gives correlations for overall force and 
moment coefficients for different angles of attack and side 
slip angles. The impact of CFD on the development of ad- 
vanced configurations is illustrated in figure 5. It shows 
the L/D performance for the configuration of figure 2 for 
across the Mach number range and compares that perfor- 
mance with existing fighters such as the F-14 and the F-15. 
A 25% to 50% increase in L/D is demonstrated. 

Figure 6 shows a complex fighter configuration with 



Fig. 3. Pressure contours at two different axial stations. 


canard, wing, vertical tail, swept-side-walled flow through 
nacelle, and a canopy. The gridding at different axial sta- 
tions is shown. Figure 7 shows the pressure contours at 
different marching planes at Moo = 2.0, a = 4°. A com- 
parison of overall force and moment coefficients is given 
in Table 2. Figure 6 illustrates the extent of geometric 
complexity that can be handled by the full potential code 
for supersonic flows. However, for transonic and subsonic 
flows where the computational domain has to extend far 
upstream and far downstream of the configuration, the re- 
quirement for a global three-dimensional grid makes treat- 
ment of complex configurations more formidable. 



Fig. 2. Computational grid for a fighter. 



x/c 


Fig. 4. Chordwise pressure distribution at 60% and 80% 
span stations; Moo = 1.6, a = 1.24°. 
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Transonic Flows — Figure 8 shows results for a canard- 
wing configuration at transonic Mach numbers. A wake 
cut is created between the trailing edge of the canard and 
the leading edge of the wing, as well as behind the trail- 
ing edge of the wing. For steady transonic computations, 
triple approximate factored time marching is performed un- 
til steady state is reached. A typical computation such as 
the one shown in figure 8 requires 100 to 200 time iterations 
requiring 60 seconds of CPU time on a CRAY-X/MP for 
80,000 grid points. 






Fig. 8. Pressure correlations for a transonic canard- 
wing configuration. 


Static Aeroelastic — Figure 9 illustrates a static aeroelastic 
computation for a flexible wing. The structural response is 
modeled using a generalized modal representation®. Within 
the aeroelastic model, a rigid wing is represented by set- 
ting the dynamic pressure to be zero. The magnitude of 
the structural deflection depends on the level of prescribed 
dynamic pressure and the generalized mode shapes. For 
an aeroelastically stable configuration, the tip load is re- 
duced once the wing undergoes static deflection. This is 
illustrated in figure 9 which shows the deflected wing shape 
along with the upper and lower surface pressures. Figure 10 
shows the Ci, versus a variation taking into account static 
flexibility. 


Moo = 1.15, a = 6“ 



Fig. 9. Static flexible computation, 



Fig. 10. Cl versus a correlation for the flexible wing. 
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Moo = 1-15 Q = 490 psf Q === 550 psf 

a - 0“ FLUTTER POINT ABOVE FLUTTER 

Q = 360 psf 
BELOW FLUTTER 

Fig. 11. Dynamic flexible computations for different dy- 
namic pressures; a) below flutter, b) at flutter 
point, and c) above flutter point. 


Dynamic Aeroelastic — Figure 11 shows dynamic flexible 
computations for three different dynamic pressure levels. 
Dynamic computations are performed as time-dependent 
calculations coupling the nonlinear aerodynamics and the 
structural response at a given time level invoking a time- 
accurate, Newton iteration procedure®. For a dynamic 
pressure level below the flutter point, the wing is aeroe- 
lastically stable as indicated by the decaying amplitude of 
oscillation in tip deflection, and tip a as a function of 
time in figure 11a. Figure 11b shows the calculation near 
flutter point with zero damping of the amplitude while fig- 
ure 11c shows results of aeroelastic divergence above the 
flutter point. 

Euler Equations 

Referring to eq. (1), the Euler equations are given by 



( ^ y 

p 


/(c-hp)u\ 

pu 

<? = 

pu 

pv 

\pw) 

,E = 

pu^ +p 
pvu 

\ pwu J 

/(c + p)v\ 


/(c + p)w\ 

pv 


pw 

puv 

,G = 

puw 

pv^ 

+p 


pvw 

V pwv / 


^ pw^ -F p > 


In the above, pressure is p computed from p = (c — p(u^ + 
+ u;^)/2))(7 — 1), density is p, Cartesian x, veloc- 
ity components are UjV,w, and the total energy per unit 
volume is e. 


In order to define the appropriate numerical fluxes 
F, G, an upwind biased scheme based on Roe’s approximate 
Riemann solver^® is employed. 

At every cell interface m-fl/2, let ^m+i/2 

denote the values of the dependent variables defined just to 
the right of and just to the left of the cell face. These 
values will be defined in the next subsection using a To- 
tal Variation Diminishing (TVD) formulation. The Rie- 
mann Solver is a mechanism to divide the flux difference 
between these neighboring states (between 

^m+ 1 / 2 ) component parts associated with each wave 
field. These can in turn be divided into those that cor- 
respond to positive and negative wave speeds. When we 
compute the numerical flux at the cell face at m + 1/2, in 
the finite-volume formulation, we will only use the cell-face 
normals defined at m -f 1/2 in the terms contribjiting to 
that representative flux. The actual fluxes E, F, G, when 
evaluated with the metrics equated to cell-face normals, can 
all be written in the same functional form given by 

F, F, G = /(Q, n,, n^, nj = /((?, N) (17) 

where the appropriate values of n;r,ny,n, are used and N 
denotes the set of those normals. Using such notation, it is 
possible to present the necessary algebra very concisely. 

Let us first denote the Jacobian matrix of the flux / 
with respect to the dependent variables Q by df/dQ. This 
Jacobian can also be called the coefficient matrix. Let us 
denote the eigenvalues of the coeflflcient matrix by A* and 
the corresponding left and right eigenvectors by £• and r*, 


9 ^ 


respectively. The matrix formed by the left eigenvectors as 
its rows is then called the left eigenvector matrix L and the 
matrix of right eigenvectors comprising the right eigenvec- 
tors as its columns is R. For our purposes, we choose an 
orthonormal set of left and right eigenvectors which implies 
that LR — RL = /, the identity matrix. In the above, the 
superscript i has been used to denote the association of the 
t-th eigenvalue with its corresponding eigenvector. Each 
eigenvalue is also associated with its own wave field. 

The underlying upwind scheme is based upon Roe’s ap- 
proximate Riemann solver. In this approach, cell interface 
values of density, velocities, and enthalpy (h = ^^/((^ - 
I)p) -h -f uP)/2) are computed using a special av- 

eraging procedure^**. 

Knowing the cell interface values, the eigenvalues and 
orthonormal set of left and right eigenvectors corresponding 
to a cell face can be computed. These may be denoted by 

^m+1/2 — ^m+l/2(Qm+l/2i ^m+l/2)i ( 18 ^) 

Kn+lj2 ~ ^m+l/2(Qm+l/2> -^m+1/2)* 

At each cell face, the positive and negative projections 
of the eigenvalues may be defined by 


where 

4 = iHQm, (Nm^l/2 + iV^^l/2)/2) (22) 

Next, we define the slope-limited values given by 

Sm+i/2 = minmod[Sj„+,/2,6S^_,/2l, 

^ <23) 

“m-i/2 = minmod[S;„_,/2,6S|„+,/2l. 

In the above, the compression parameter 6 is to be taken as 
the following function of the accuracy parameter <p which 
is explained shortly. 


The minmod slope-limiter operator is 

minmod[x, j/| = sign(x) max[0,min{|x|, j/sign(x)}] (25) 

Then, the left state at the cell interface at m 4- 1/2 and 
the right state at the cell interface m — 1/2 can be defined 
to be 


i ^ 

14 -^^,- ^ 1 — ^ 

4 «m+l/2+ 4 ®m-l/2^ 

)rU 

^m-1/2 “ “ Xrf f 

1 + ^ 1 — ^ ^ 
4 «m-l/2 + 4 

)rU 

t 


(26) 

where 



rL = rHQ.n. 

1 (^m+1/2 + ■Wm-l/2)/2) 

(27) 


At maxima and minima, the minmod operator returns a 
zero value and the left and right states reduce to 


. (186) 

Now, the numerical flux /m+1/2 is constructed from 
/m+1/2 

= 2 [*^(^m+l/2>^Wl/2) +/(Qm+l/2>^m+l/2)] 
^m+1/2 ~ '^m+1/2) ^m+l/2^m+l/2 

~ <^(^m-M/ 2 > ^rn+112) + ^ ^ ‘^m+l/ 2 ^m-|-l/ 2 ^m-|-l /2 

% 

(19) 

In the above equation, 

<^m+l /2 = ^+l/ 2 (Qm+l /2 ^m+ 1 / 2 )* 

We can construct upwind-biased schemes of varying 
accuracies by properly defining the left and right states 
used in the last subsection. We present here a family of 
schemes. For use in what follows, let us now define 

^m+l/2 ~ ~ Qm)i 


= Qm 

which result in a first-order accurate scheme locally. 

More details on this Euler solver can be found in 
Refs. 12-14. Now some results are presented to illustrate 
the unified Euler solver capability. 

Euler Results 

Supersonic Flows — Figure 12a shows an elliptic waverider 
geometry typical of hypersonic configurations. Typically, 
waveriders are designed to have a lift-producing lower sur- 
face with a freestream aligned upper surface. Figure 12b 
shows Mach number contours at different Mach numbers 
and angles of attack. At the design point (Mo© = 4, a = 
0®), the shock is at the leading edge while at off-design flow 
conditions the shock moves away from the leading edge. 
The upwind, TVD based Euler solver implemented in the 
code does not exhibit any numerical instability problems in 
capturing strong shocks. Figures 12c and 12d show compar- 
isons of surface pressures and pitching moment coefficients 
with experimental data^^ and other available methods^^. 
The full potential method compares well with the Euler re- 
sults when the shock is weak. The shock strength starts to 
become more pronounced for M^o > 4, a > 5® as indicated 
by the deviation of the isentropic full potential results from 
the correct Euler solutions. 
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Fig. 12. Euler results for a waverider configuration; a) ge- 
ometry and gridding, b) Mach contours at differ- 
ent Mach numbers, c) lower surface pressures, 
and d) pitching moment correlation. 


Multizone Computation — For treatment of complex three- 
dimensional multibody flows, the gasdynamic solvers are 
provided with a multizonal capability where the physical 
domain of interest is subdivided into multizones requiring 
single gridding procedures within each zone. Across the 
zonal boundaries proper flux balancing is maintained to 
avoid spurious numerical errors originating at the interface. 
The zonal interface can be permeable or impermeable and 
can also be a boundary of flow discontinuity such as a shock 
or sonic surface. 

Figure 13 shows the Space Shuttle mated configuration 
with the Orbiter mounted on top of the External Tank and 
the Solid Rocket Boosters. A single zone gridding that 
can treat every component of this multibody as a constant 


coordinate surface, though possible, can be cumbersome 
to construct. A five-zone gridding in the axial plane is 
generated to study this multibody problem at supersonic 
Mach numbers. Pressure contours and gridding are shown 
at different marching stations for Moo = 1*8, Of = 0°. The 
presence of a shock around the Shuttle OMS pod (station 
C) is clear. 

TV*ansonfc Flow — Figure 14 shows transonic results for the 
ONERA-M6 wing. The double shock pattern on the upper 
surface at Moo “ 0.84, a = 3.06° is well captured by the 
Euler code. 
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MULTIZONE COMPUTATION 

Moo = 1.8, a = 0“ 



Fig. 13. Multizone treatment of the Shuttle mated con- 
figuration; Moo = 1-8, Of = 0®. 


Navier- Stokes Equations 

The purpose in developing powerful, robust and effi- 
cient Euler solvers is not just to study inviscid strongly 
shocked, rotational flows, but also to use them as a step- 
ping stone in devising Navier-Stokes methods for solving 
viscous flow problems. Many problems of real interest in 
advanced aerospace and configuration development do re- 
quire the use of Navier-Stokes methods. Some of the flows 
that can only be modeled using the Navier-Stokes equations 
are: 

1) attached flows with a) tip vortex, b) wing/body junc- 
tion vortex, and c) cross-flow leading edge vortex; 

2) separated flows (leading edge separation and shock- 
boundary layer separation); 

3) acoustics/unsteady phenomena (cavity flow and inter- 
nal flow-induced vibrations); 

4) high Mach number flows with significant heating; and 

5) reacting flows (combustion involving chemical kinet- 
ics). 

Referring to eq. (1), the Reynolds- averaged form of the 
Navier-Stokes equations is represented by 



Fig. 14. Transonic results for the ONERA-M6 wing; 
Moo =0.84, a = 3.06®. 
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(29) 

where E,*n, and Gin are the inviscid Euler fluxes given 
by eq. (16) and the viscous fluxes E„, F„, and are given 
by 
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Here, Re is the Reynolds number, Pr is the Prandtl num- 
ber, K is the thermal conductivity, and T is the specific 
internal energy given by T = p/[(Tf — l)p]- The terms r,*, 
T*y, »•**, »•»*- »•»». ^vz, fzx, r^y, and t„ are given by 
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(31) 


where 


5u dv dw. 
dx dy dz 


and fi is the coefficient of viscosity. 

The viscosity coefficient for turbulent flows is modeled 
as the sum of the laminar and turbulent viscosities in the 
eddy viscosity approach. The turbulent eddy viscosity is 
usually computed using one of two popular techniques, 1) 
by using the Baldwin-Lomax or other algebraic eddy vis- 
cosity formulation, and 2) by using a two-equation model 
such as the A; - € formulation. 


The k-e model often used is the standard high Reynolds 
number form of the equations. Even though the k-e model 
can take more time to solve than the simpler algebraic eddy 
viscosity models this is justifiable since the k-e model is gen- 
erally applicable to a much wider class of flows. The kinetic 
energy equation is derived from the Navier-Stokes equations 
with the main limiting criterion being that it assumes local 
isotropy. The dissipation equation is not exact but is mod- 
elled to represent physical processes similar to those of the 
kinetic energy equation. Even with these assumptions the 
k-e equations have a proven capability of adequately pre- 
dicting a large range of complex flows, including anisotropic 
ones. 


The k-e equations may be solved using the same up- 
wind, TVD formulations applied to the Euler and Navier- 
Stokes equations. Referring to eq. (1), the k-e equations 
can be written as 


-it)’ 
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pvk 
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~ Re dy / V Re dz / 
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The only exception is that the k-e equations, in addition 
to the above, involve a source term on the right hand side 
given by 

' P-peRe \ 


In 


^ Re 


and f32l. 
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(33) 


fXk = + 

/X, = (/X-h/itM) 

k is kinetic energy, e is turbulent dissipation, and fit is 
turbulent eddy viscosity. P represents the production of 
kinetic energy and the following simplified form of it is used 

P = fit(ul + Wy d- ul). (34) 

The k-e model still employs the eddy viscosity /difTusiv- 
ity concept as it relates eddy viscosity to the kinetic energy 
and dissipation by 

fit = C^ip—. (35) 

This eddy viscosity is then used to create an effective vis- 
cosity (m + Mt) which replaces fi in the Reynolds- averaged 
Navier-Stokes equations. To solve the above turbulence 
model the following constants must be specified: <Jk = EO, 
(T^ = 1.3, Cl = 1.44, C 2 = 1.92, and C^ 0.09, 

One of the highlights of the Navier-Stokes activity is 
the modeling of turbulence for separated flows. The new 
turbulence model is based on experimental observations of 
separated turbulent flows. The model prescribes turbulence 
kinetic energy (k) and its dissipation (e) analytically inside 
separation bubbles. A Gaussian variation of k normal to 
walls is assumed. The length scale of turbulence within 
bubbles is proportional to the local distance from the wall to 
the edge of the viscous sublayer, which is located outside the 
backflow region, as shown in figure 15. The latter feature 
is a basic assumption of the model. 

The stress scale is the local maximum Reynolds stress, 
which typically occurs around the middle of the boundary 
layer, well outside the bubble. This scale must be supplied 
by the turbulence model used beyond separated regions. 

The main equations of the model are given in Ref. 23. 
A simple formula for eddy viscosity distribution within the 
separation bubble results, and is used to provide eddy vis- 
cosity for the Reynolds-averaged equations when perform- 
ing the calculations inside the bubble. Outside of it, an- 
other turbulence model (e.g., Baldwin-Lomax or Ar-e) sup- 
plies the values of eddy viscosity. 
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Fig. 15. Schematic view of separated flow bubble and ba- 
sic nomenclature. 


Navier-Stokes Results 

The algebraic k~e turbulence model for separated 
flows, in conjunction with the Baldwin-Lomax model, has 
been incorporated into a finite volume, time-marching, 
multizonal Navier-Stokes code^^^^®, featuring an implicit 
upwind-biased scheme, approximate factorization, and To- 
tal Variation Diminishing discretization for high accuracy. 
When a separation bubble exists, the standard Baldwin- 
Lomax or the k-e model is used to compute eddy viscos- 
ity outside the backflow region, while the separation model 
(figure 15) provided eddy viscosity within the separation 
bubble. 

Several unit problems are computed to check the va- 
lidity of the Navier-Stokes code with a turbulence model 
for treating separated flows. 

As a first computational test of the new turbulence 
model, a transonic flow calculation over an axisymmetric 
boattail with a cylindrical extension (solid plume simulator) 
has been performed. Figure 16a shows the chosen geome- 
try. This case involves a moderate-sized separation bubble 
at the end of the boattail. The data of Ref. 32 at Afo© = 0.8 
and a Reynolds number of 1.8 x 10®, based on maximum 
model diameter, were used for comparisons with the calcu- 
lations. A 65 X 40 grid was employed, with 23 points normal 
to the wall lying inside the separation bubble at the loca- 
tion of its maximum height. Figure 16b shows a detail of 
the computational mesh. 

In figure 16c, pressure coefficient at the wall, calcu- 
lated using the new separation turbulence model, is com- 
pared with experimental data of Ref. 32 as well as with a 
calculation which used the Baldwin-Lomax model by itself. 
The figure also indicates the location and extent of the sep- 
arated region. The advantage of using the separation model 
is demonstrated by the significant improvement in predict- 
ing the pressure through the separated zone, compared to 
the corresponding calculation without the model. 

Figure 16d compares skin friction distribution, as cal- 
culated using the new model, with the corresponding cal- 
culation done without it. A larger separation bubble is 
predicted by the former. No data are available for compar- 
ison. 
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The second test of the separation model was chosen to 
be a supersonic flow over a 24° compression ramp. Detailed 
experimental data are available for this case^^“®^, which 
involves a rather large separated flow region. The 65 x 
35 grid used for the calculations is shown in figure 17a. 
About 75% of the normal-to- wall mesh points were located 
within the shear layer. This case was run at Moo = 2.85 
and a Reynolds number based on incoming boundary layer 
thickness of 1.6 x 10®. 

Figure 17b shows a comparison between calculations 
and experimental data of wall pressure distribution. The 
improved predictional capability due to inclusion of the sep- 
aration model is evident. 

Figure 17c compares calculations with experimental 
data of skin friction. The advantage of using the sepa- 
ration model over the regular Baldwin-Lomax model is ev- 
ident throughout the separation bubble. Incorporation of 
the new model enables precise prediction of reattachment 
location, a difficult task for this flow. 

A streamline plot, resulting from the calculation with 
the separation model, is shown in figure 17d. The predicted 
extent of the separated region agrees quite well with the 
experimentally observed locations of separation and reat- 
tachment, also indicated in the figure. 

The third test of the new backflow model was the 
backward-facing step case reported in Ref. 36, with an 
inflow Mach number of 0.128 and a Reynolds number of 
31,250 based on inflow conditions and the step height as a 
reference length. 

Figure 18a shows the geometry and the two-zone com- 
putational grid, with a 42 x 22 mesh used in the subdomain 
above the step, and a 36 x 20 mesh used in the subdomain 
downstream of the step. 

Pressure distribution along the step-side wall is shown 
in flgure 18b, as resulting from the current approach and 
from Sindir’s^^ calculations using the k-e model. Compar- 
ison with the data indicates a slight advantage in using the 
new algebraic model over the k-e model. 

Figure 18c shows skin friction distribution on the step- 
side wall. The new backflow model enables improved pre- 
diction and significantly better performance in the reat- 
tachment zone. In the vicinity of the step corner, the skin 
friction is positive, indicating a small counter-rotating vor- 
tex, in agreement with the data. 

In figure 18d, streamwise velocity profiles at two loca- 
tions are shown, one upstream of the reattachment region, 
the other downstream of it. Agreement with the data is 
very good at the former location, where the flow is sepa- 
rated, and fair at the latter, where a somewhat sluggish 
boundary layer recovery is predicted for the lower part of 
the profile. 

Figure 18e shows Reynolds stress profiles at the loca- 
tions corresponding to those of figure 18d. While the shape 
of the calculated profile agrees with the experimental one 
at the upstream location, the magnitude is overpredicted 
roughly by a factor of two. In the downstream location, 
however, agreement with data is quite good, although the 
lower part of the calculated profile is again overpredicted. 
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Since store separation and aerodynamic drag due to 
open cavities are important issues in design, cavity compu- 
tations are vital for validation and prediction. A laminar 
three-dimensional cavity computation has been done. The 


specific cavity is a simplified version of the F-111 weapons 
bay. Figure 19a shows velocity directions down the center- 
line of the cavity. Figure 19b shows the velocity directions 
of the secondary motion on a cross plane of the cavity. More 
cavity results can be found in Ref. 38. 
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Fig. 19. F-111 cavity flow; a) center plane, b) cross plane. 
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Fig. 20c. Pressure contours 



Fig. 20d. Wall pressure comparison 



Fig. 20e. Skin friction 



Fig. 20f. Plume development 


Figure ?0a shows an axisymmetric nozzle and the cor- 
responding two-zone gridding is given in figure 20b. The 
external inflow was fixed at Moo — 20, while the nozzle in- 
flow Mach number was varied from 4.0 to 1.5. The nozzle 
inflow pressure was 217 times that of the external inflow. 
For the viscous calculations, the Reynolds number based 
on nozzle inflow conditions was varied from 6 x 10^ /ft at 
Mach 4 to 2.25 x 10^/ft at Mach 1.5. All solid surfaces were 
held at a constant wall temperature Tw = 0.34Too,noz- Cal- 
culations for this case were run using both the Reynolds- 
averaged Navier-Stokes code (RANS) and the Euler code. 
Figure 20c shows pressure contours for the Navier-Stokes 
computation. Figure 20d compares pressure distributions 
on the nozzle wall and on the outer plate surface plus the 
wake region downstream of it, as resulting from the RANS 
and Euler calculations. Except for the milder corner ex- 
pansion due to viscous displacement effects, the two calcu- 
lations predict the same nozzle wall pressure distribution. 
In the near wake region downstream of the splitter plate, 
however, the interaction between the shear layers from both 
sides of the plate and the plume-induced shock/boundary 
layer interaction in the vicinity of the plate trailing edge 
modify the pressure distribution as compared with the in- 
viscid prediction. Further downstream the two predictions 
coincide. Figure 20e shows skin friction on the two walls, in- 
dicating no separation of the boundary layers. The stream- 
line plot in figure 20f shows the plume development as pre- 
dicted by the RANS code. 

Figure 21a shows the geometry and the multizone grid- 
ding for a ramjet configuration. Figure 21b shows viscous 
Mach number contours for an inflow Mach number of 4.03. 

The Navier-Stokes code developed at the Rockwell Sci- 
ence Center is still undergoing validation tests on several 
unit problems for possible improvements in modeling tur- 
bulence of separated flows. Future applications will involve 
wings, wing- body combinations at high a, and cavity-store 
acoustics and separation studies. 
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Laser-Material Interaction 

The paper so far dealt with the development and ap- 
plication of computational algorithms for solving aerody- 
namic problems. Extension of these CFD methods to solv- 
ing problems in other disciplines that are governed by an 
appropriate set of partial differential equations is very at- 
tractive. One such application is to study the problem of 
heat transfer in materials subjected to intense laser heating. 

Laser heat treatment of materials (especially iron-base 
alloys and carbon-carbon composites) for various industrial 
applications is becoming very attractive due to ease in the 
controllability and in generation of laser beams. For ex- 
ample, use of laser as a heat source in enhancing materi- 
als resistance to surface wear and corrosion through solid 
state phase transformations (without melting) and rapid 
solidification (with shallow melting), in achieving a desired 
homogeneous molten weld pool, and in obtaining a unique 
surface composition through coating or cladding, as a vi- 
able economical process, has been well proven in laboratory 
settings. Transitioning this process technology from a labo- 
ratory setting to an industrial environment requires a better 
understanding of the role of various controlling parameters, 
such as the cross-section of the laser, power intensity of the 
laser, velocity of the moving laser or the workpiece, and the 
material properties themselves in determining the quality 
of the surface modification process. Optimization of these 
controlling process parameters through theoretical model- 
ing and computational simulation can lead to achieving the 
desired properties of surface treatment. 

Figure 22 shows the schematic of a laser melted ma- 
terial pool. When the workpiece is swept under the beam, 
a self quenched heat treated zone is obtained along the 
surface. Dimensions of the melted zone (T > 7V„, where 
Tm is the melting temperature) and the heat affected zone 
(T > Tc) are controlled by absorbed laser beam power den- 
sity, beam size, and travel speed. 



MZ - MELTED ZONE 

HAZ - HEAT AFFECTED ZONE 

Tjj - PHASE TRANSFORMATION TEMP. 

T^ - MELTING TEMP. 
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Fig. 22. Surface tension induced convective heat transfer. 


The high surface temperature gradients in the melted 
zone create a variation of surface tension along the surface 
which is balanced by shear forces. This balancing shear 
force is created by setting up a counter rotating vortex flow 
within the molten zone as described in figure 22. 

The physics of modeling the heat transfer process oc- 
curring in figure 22 involves both conduction and convec- 
tion in a high gradient thermal field. The equation that 
best describes the physical phenomena of this problem is 
the incompressible Navier-Stokes equations. Referring to 
eq. (3), the following describes the coupling between the 
temperature field and the velocity field. 
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a = 


where r,y = u Modeling of turbulence in r,-j 

is neglected in the present formulation and only the laminar 
stress tensor is considered. 


In these equations, Cp is the specific heat, p is the den- 
sity, k is the thermal conductivity, T is the temperature, 
and j/ is the kinematic viscosity. The induced velocity field 
in the molten pool is represented by uy. 

In the unmelted region of the material where the heat 
transfer process is purely due to conduction (temperatures 
below melting), only the energy equation needs to be solved 
for temperature (see Ref. 39). 

On the outer surface of the molten pool, the force bal- 
ance equations are 

du ,dT 

dl ,dT (37) 

^dz dy 

where g* is the rate of change of surface tension (TV/m/A:) 
with temperature and p, is the coefficient of viscosity of the 
molten pool (N sec/m^; N is Newton). 

At the liquid-solid interface u = u = ty = 0 and T = 
Tfn- More details on the boundary condition can be found 
in Ref. 40. 


The computational method employs an implicit triple 
approximate factorization scheme to solve the energy equa- 
tion in terms of temperature and an explicit treatment for 
the three momentum equations and the continuity equa- 
tion. The pressure field is updated at each time level using 
a Poisson solver to satisfy the continuity equation. 
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Laser Results 

A sample result for a rectangular workpiece undergoing 
melting is presented. 

Figure 23 shows results for a typical case (5000 W laser, 
4x3 mm laser cross section, 1.26 mm/sec beam travel 
speed). A point in the workpiece is considered heat af- 
fected if that point experienced temperatures above 750° C 
and the melt zone corresponds to temperatures > 1500°C. 
At a given instance of time (laser beam around the hallway 
length of the workpiece), the instantaneous temperature 
distribution, melt zone shape, and the three-dimensional 
vortex mixing material flow are shown in the figure. The 
cross-sectional view, longitudinal view (plane of symme- 
try), and the top view of the melted zone clearly reveal the 
convection process induced by the surface tension driven 
flow. Pure conduction treatment of this heat transfer prob- 
lem (no convection model, i,e,, surface tension gradients 
set to zero) results in temperature levels not comparable 
with experimental observations. The computational proce- 
dure of the present study incorporating the surface tension 
driven convective heat transfer process produces melt zone 
and heat affected zone shapes very similar to experimental 
data. More results are presented in Ref. 40. 

Application of this work to study problems of deep 
welding (solute redistribution, microstructure of the heat 
affected zone, and residual stress state) and characteriza- 
tion of the surface ripples are some of the ongoing projects. 


Electromagnetic Scattering 

The objective is to develop time-dependent finite dif- 
ference methods to solve the Maxwell equations to study 
the problem of electromagnetic scattering from dielectric 
and perfectly reflecting objects. Although techniques based 
on the integral form of the equations are available, they are 
usually restricted in their application due to various simpli- 
fications made in the formulation. Solution techniques to 
the differential equation usually provide a general-purpose 
capability with fewer restrictions than techniques based on 
the integral approach. Based on proven CFD methods, it is 
desirable to develop an efficient finite difference technique 
for the Maxwell equations. 

Referring to eq. (1), the Maxwell equations take the 




Fig. 23. Recirculating flow within the molten pool; beam 
power = 5000 W, beam shape = 4x3 mm, pro- 
cess speed = 1.27 mm/sec. 



where c^, Cy, and are the electric field components along 
X, y, and z, and similarly Hx, Hy, and Hz are the mag- 
netic field components. The parameters e, /i, and a rep- 
resent permittivity, permeability, and conductivity of the 
medium through which the electromagnetic wave is propa- 
gating. Equation (38) is hyperbolic and has real eigenval- 
ues and a linearly independent set of eigenvectors. Upwind 
schemes developed for the Euler equations are ideal for solv- 
ing the Maxwell equations. The objective is to solve eq. (3) 
subject to an incident wave to compute the equivalent sur- 
face current on the object given hy n y. H where n is the 
surface normal and H is the magnetic field vector. Once 
the equivalent surface current is known on the dielectric 
object and around any contour encompassing the object, 
the radar cross section (RCS) information can be obtained 
using a near field-to-far field transformation (Refs. 41,42). 
The RCS information depends on the intensity of the scat- 
tered wave. 

RCS Results 

The Maxwell equations in two dimensions can be spe- 
cialized for a transverse magnetic (TM) wave (e* = Cy = 
= 0), or for a transverse electric (TE) wave (Hx - 
Hy = 0,c, =0). Development of computational algorithms 
for studying electromagnetic scattering from perfectly con- 
ducting or dielectric objects can benefit from solving the 
TM or the TE wave problem. Preliminary results are re- 
ported here for the electromagnetic scattering from a per- 
fectly conducting square cylinder for an incident plane wave 
of the form c,- = Co5»n^(rco5(^ - ~ c<) where Cq is 

the amplitude, k = 2ir/A, A is the wavelength, c is the 
wave speed, is the angle of incident wave with respect 
to the I- axis. Figures 24a and 24b show the surface cur- 
rent Jz — n X H on the cylinder for two different incident 
wave angles. The correlation of n x R obtained using a 
simple, first-order accurate, upwind, explicit scheme with 
an existing method known as method of moments (MOM) 
is good. Knowing this J, information, the RCS value for 
different viewing angles can be computed. Development 
of higher order accurate, upwind schemes based on Euler 
solvers is currently in progress. Some of the numerical is- 
sues to be addressed in this development are 1) higher order 
accurate nonreflecting farfield conditions based on charac- 
teristic theory, 2) grid resolution requirements for high fre- 
quency (small A or large k) incident waves, 3) boundary 
condition treatment for radar absorbing materials taking 
into account frequency dependence on €, /i and <r, 4) multi- 
zone gridding techniques for interior and external regions, 
and 5) near-field to far-field transformations to derive RCS 
values from the near-field nx H and nx E information. 
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Fig. 24a Surface current on a square cylinder 
for a right moving incident wave. 



Fig. 24b Surface current on a square cylinder 
for an incident wave at 45® . 



CONCLUSIONS 

The state of the art of Computational Fluid Dynamics 
has taken rapid strides in recent years with the development 
and application of unified, robust, and efficient methods. 
The advances in CFD are also beginning to make a positive 
impact on other areas of mathematical science, leading to 
the emergence of the concept of “Computational Science”. 
In this new spirit, this paper has presented a unification of 
algorithms and their application to fluid dynamics, electro- 
magnetics, and material characterization. 
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